/*
* ManagedLapack.cs
*
* Managed provider for LAPACK routines.
*
* Many of the routines are a port of JAMA and Jampack code.
*
* Copyright (c) 2005, dnAnalytics. All rights reserved.
*/
using dnAnalytics.LinearAlgebra.Decomposition;

namespace dnAnalytics.LinearAlgebra
{
    internal sealed class ManagedLapack : ISimplifiedLapack
    {
        #region fields
        private ISimplifiedBlas blas = new ManagedBlas();
        #endregion fields
        #region Public Methods
        public void LUFactor(int order, double[] factor, int[] pivots)
        {
            for (int i = 0; i < order; i++)
            {
                pivots[i] = i;
            }
            //double[] LUrowi = new double[order];
            double[] LUcolj = new double[order];

            int indexj, indexjj, kmax, indexk, indexkp, indexkj;
            double temp;
            // Outer loop.
            for (int j = 0; j < order; j++)
            {
                indexj = j * order;
                indexjj = indexj + j;
                // Make a copy of the j-th column to localize references.
                for (int i = 0; i < order; i++)
                {
                    LUcolj[i] = factor[indexj + i];
                }

                // Apply previous transformations.
                for (int i = 0; i < order; i++)
                {
                    // Most of the time is spent in the following dot product.
                    kmax = System.Math.Min(i, j);
                    double s = 0.0;
                    for (int k = 0; k < kmax; k++)
                    {
                        s += factor[k * order + i] * LUcolj[k];
                    }

                    factor[indexj + i] = LUcolj[i] -= s;
                }

                // Find pivot and exchange if necessary.
                int p = j;
                for (int i = j + 1; i < order; i++)
                {
                    if (System.Math.Abs(LUcolj[i]) > System.Math.Abs(LUcolj[p]))
                    {
                        p = i;
                    }
                }
                if (p != j)
                {
                    for (int k = 0; k < order; k++)
                    {
                        indexk = k * order;
                        indexkp = indexk + p;
                        indexkj = indexk + j;
                        temp = factor[indexkp];
                        factor[indexkp] = factor[indexkj];
                        factor[indexkj] = temp;
                    }
                    pivots[j] = p;
                }

                // Compute multipliers.

                if (j < order & factor[indexjj] != 0.0)
                {
                    for (int i = j + 1; i < order; i++)
                    {
                        factor[indexj + i] /= factor[indexjj];
                    }
                }
            }
        }

        public void LUInverse(int order, double[] factor, int[] pivots)
        {
            double[] ident = new double[order * order];
            for (int i = 0; i < order; i++)
            {
                ident[i * order + i] = 1;
            }

            LUSolve(order, order, factor, pivots, ident);

            System.Buffer.BlockCopy(ident, 0, factor, 0, ident.Length * Constants.SizeOfDouble);
        }

        private static double[] Pivot(int m, int n, double[] B, int[] pivots)
        {
            double[] ret = new double[m * n];
            System.Buffer.BlockCopy(B, 0, ret, 0, B.Length * Constants.SizeOfDouble);

            int indexk, indexkp, indexkj;
            double temp;

            for (int i = 0; i < pivots.Length; i++)
            {
                if (pivots[i] != i)
                {
                    int p = pivots[i];
                    for (int j = 0; j < n; j++)
                    {
                        indexk = j * m;
                        indexkp = indexk + p;
                        indexkj = indexk + i;
                        temp = ret[indexkp];
                        ret[indexkp] = ret[indexkj];
                        ret[indexkj] = temp;
                    }
                }
            }
            return ret;
        }

        public void LUSolve(int order, int columns, double[] factor, int[] pivots, double[] rightSide)
        {
            double[] data = Pivot(order, columns, rightSide, pivots);

            // Solve L*Y = B(piv,:)
            int korder = 0;
            int index = 0;
            for (int k = 0; k < order; k++)
            {
                korder = k * order;
                for (int i = k + 1; i < order; i++)
                {
                    for (int j = 0; j < columns; j++)
                    {
                        index = j * order;
                        data[i + index] -= data[k + index] * factor[i + korder];
                    }
                }
            }
            // Solve U*X = Y;
            for (int k = order - 1; k >= 0; k--)
            {
                korder = k + k * order;
                for (int j = 0; j < columns; j++)
                {
                    data[k + j * order] /= factor[korder];
                }
                korder = k * order;
                for (int i = 0; i < k; i++)
                {
                    for (int j = 0; j < columns; j++)
                    {
                        index = j * order;
                        data[i + index] -= data[k + index] * factor[i + korder];
                    }
                }
            }
            System.Buffer.BlockCopy(data, 0, rightSide, 0, data.Length * Constants.SizeOfDouble);
        }

        public double Norm(char norm, int rows, int columns, double[] data)
        {
            double ret = 0;
            switch (norm)
            {
                case '1':
                case 'O':
                case 'o':
                    for (int j = 0; j < columns; j++)
                    {
                        double s = 0;
                        for (int i = 0; i < rows; i++)
                        {
                            s += System.Math.Abs(data[j * rows + i]);
                        }
                        ret = System.Math.Max(ret, s);
                    }
                    break;

                case 'I':
                case 'i':
                    {
                        for (int i = 0; i < rows; i++)
                        {
                            double s = 0;
                            for (int j = 0; j < columns; j++)
                            {
                                s += System.Math.Abs(data[j * rows + i]);
                            }
                            ret = System.Math.Max(ret, s);
                        }
                    }
                    break;

                case 'F':
                case 'f':
                case 'E':
                case 'e':
                    {
                        double[] trans = new double[data.Length];
                        for (int i = 0; i < rows; i++)
                        {
                            for (int j = 0; j < columns; j++)
                            {
                                trans[i * columns + j] = data[j * rows + i];
                            }
                        }

                        double[] tmp = new double[rows * rows];
                        blas.Multiply(rows, rows, columns, data, trans, tmp);
                        for (int i = 0; i < rows; i++)
                        {
                            ret += System.Math.Abs(tmp[i * rows + i]);
                        }
                        ret = System.Math.Sqrt(ret);
                    }
                    break;

                default:
                    throw new System.NotSupportedException();
            }
            return ret;
        }

        public double ConditionNumber(char norm, int order, double anorm, double[] luData, int[] pivots)
        {
            double[] inverse = new double[luData.Length];
            System.Buffer.BlockCopy(luData, 0, inverse, 0, luData.Length * Constants.SizeOfDouble);
            LUInverse(order, inverse, pivots);

            double invNorm = 0;
            switch (norm)
            {
                case '1':
                case 'O':
                    invNorm = Norm('1', order, order, inverse);
                    break;
                case 'I':
                    invNorm = Norm('I', order, order, inverse);
                    break;

                default:
                    throw new System.NotSupportedException();
            }
            return 1.0 / (anorm * invNorm);
        }

        public int CholeskyFactor(int order, double[] data)
        {
            int i, j, k, tmp, index;
            double s, d;
            double[] factor = new double[data.Length];
            for (j = 0; j < order; j++)
            {
                d = 0.0;
                for (k = 0; k < j; k++)
                {
                    s = 0.0;
                    for (i = 0; i < k; i++)
                    {
                        s += factor[i * order + k] * factor[i * order + j];
                    }
                    tmp = k * order;
                    index = tmp + j;
                    factor[index] = s = (data[index] - s) / factor[tmp + k];
                    d += s * s;
                }
                index = j * order + j;
                d = data[index] - d;
                if (d <= 0.0)
                {
                    return 1;
                }
                factor[index] = System.Math.Sqrt(d);
                for (k = j + 1; k < order; k++)
                {
                    factor[k * order + j] = 0.0;
                }
            }
            System.Buffer.BlockCopy(factor, 0, data, 0, factor.Length * Constants.SizeOfDouble);
            return 0;
        }

        public void CholeskySolve(int order, int columns, double[] data, double[] rhs)
        {
            double sum = 0;
            int cindex = 0;
            int iindex = 0;
            for (int c = 0; c < columns; c++)
            {
                cindex = c * order;
                // Solve L*Y = B;
                for (int i = 0; i < order; i++)
                {
                    sum = rhs[c * order + i];
                    for (int k = i - 1; k >= 0; k--)
                    {
                        sum -= data[k * order + i] * rhs[cindex + k];
                    }
                    rhs[cindex + i] = sum / data[i * order + i];
                }

                // Solve L'*X = Y;
                for (int i = order - 1; i >= 0; i--)
                {
                    sum = rhs[cindex + i];
                    iindex = i * order;
                    for (int k = i + 1; k < order; k++)
                    {
                        sum -= data[iindex + k] * rhs[cindex + k];
                    }
                    rhs[cindex + i] = sum / data[iindex + i];
                }
            }
        }

        public void QRFactor(int m, int n, double[] q, double[] r, double[] tau)
        {
            for (int i = 0; i < m; i++)
            {
                q[i * m + i] = 1;
            }
            int minmn = m < n ? m : n;

            double[][] u = new double[minmn][];
            for (int i = 0; i < minmn; i++)
            {
                u[i] = Householder.GenerateColumn(r, m, i, m - 1, i);
                Householder.UA(u[i], r, m, i, m - 1, i + 1, n - 1);
            }

            for (int i = minmn - 1; i >= 0; i--)
            {
                Householder.UA(u[i], q, m, i, m - 1, i, m - 1);
            }
        }

        public void QRSolve(int m, int n, int bn, double[] q, double[] r, double[] b, double[] tau)
        {
            // Compute Y = transpose(Q)*B
            double[] column = new double[m];
            int jm, im, km;
            double s;
            for (int j = 0; j < bn; j++)
            {
                jm = j * m;
                for (int k = 0; k < m; k++)
                {
                    column[k] = b[jm + k];
                }
                for (int i = 0; i < m; i++)
                {
                    s = 0;
                    im = i * m;
                    for (int k = 0; k < m; k++)
                    {
                        s += q[im + k] * column[k];
                    }
                    b[jm + i] = s;
                }
            }

            // Solve R*X = Y;
            for (int k = n - 1; k >= 0; k--)
            {
                km = k * m;
                for (int j = 0; j < bn; j++)
                {
                    b[j * m + k] /= r[km + k];
                }
                for (int i = 0; i < k; i++)
                {
                    for (int j = 0; j < bn; j++)
                    {
                        jm = j * m;
                        b[jm + i] -= b[jm + k] * r[km + i];
                    }
                }
            }
        }

        //port of th linpack SVD routing, we should switch to the lapack SVD routine
        [System.Diagnostics.CodeAnalysis.SuppressMessage("Microsoft.Maintainability", "CA1502:AvoidExcessiveComplexity")]
        public int SVDDecomposition(int rows, int cols, double[] matrix, double[] s, double[] u, double[] v)
        {
            int MAXITER = 1000;
            // Derived from LINPACK code.
            // Initialize.
            int mn = System.Math.Min(rows, cols);
            bool computeVectors = true;
            if (u == null)
            {
                computeVectors = false;
            }
            double[] e = new double[cols];
            double[] work = new double[rows];

            int i, iter, j, k, kase, l, lp1, ls = 0, lu, m, nct, nctp1, ncu, nrt, nrtp1;
            double b, c, cs = 0.0, el, emm1, f, g, scale, shift, sl,
                sm, sn = 0.0, smm1, t1, test, ztest, xnorm, enorm;
            double t, r;

            ncu = rows;

            //   reduce matrix to bidiagonal form, storing the diagonal elements
            //   in s and the super-diagonal elements in e.
            int info = 0;
            nct = System.Math.Min(rows - 1, cols);
            nrt = System.Math.Max(0, System.Math.Min(cols - 2, rows));
            lu = System.Math.Max(nct, nrt);

            for (l = 0; l < lu; l++)
            {
                lp1 = l + 1;
                if (l < nct)
                {
                    // compute the transformation for the l-th column and
                    // place the l-th diagonal in s[l].
                    xnorm = Dnrm2Column(matrix, rows, l, l);
                    s[l] = xnorm;
                    if (s[l] != 0.0)
                    {
                        if (matrix[l * rows + l] != 0.0)
                        {
                            s[l] = Dsign(s[l], matrix[l * rows + l]);
                        }
                        DscalColumn(matrix, rows, l, l, 1.0 / s[l]);
                        matrix[l * rows + l] = 1.0 + matrix[l * rows + l];
                    }

                    s[l] = -s[l];
                }

                for (j = lp1; j < cols; j++)
                {
                    if (l < nct)
                    {
                        if (s[l] != 0.0)
                        {
                            // apply the transformation.
                            t = -Ddot(matrix, rows, l, j, l) / matrix[l * rows + l];
                            for (int ii = l; ii < rows; ii++)
                            {
                                matrix[j * rows + ii] += t * matrix[l * rows + ii];
                            }
                        }
                    }

                    //place the l-th row of matrix into  e for the
                    //subsequent calculation of the row transformation.
                    e[j] = matrix[j * rows + l];
                }

                if (computeVectors && l < nct)
                {
                    // place the transformation in u for subsequent back multiplication.
                    for (i = l; i < rows; i++)
                    {
                        u[l * rows + i] = matrix[l * rows + i];
                    }
                }

                if (l < nrt)
                {
                    // compute the l-th row transformation and place the l-th super-diagonal in e(l).
                    enorm = Dnrm2Vector(e, lp1);
                    e[l] = enorm;
                    if (e[l] != 0.0)
                    {
                        if (e[lp1] != 0.0)
                        {
                            e[l] = Dsign(e[l], e[lp1]);
                        }
                        DscalVector(e, lp1, 1.0 / e[l]);
                        e[lp1] = 1.0 + e[lp1];
                    }
                    e[l] = -e[l];

                    if (lp1 < rows && e[l] != 0.0)
                    {
                        // apply the transformation.
                        for (i = lp1; i < rows; i++)
                        {
                            work[i] = 0.0;
                        }

                        for (j = lp1; j < cols; j++)
                        {
                            for (int ii = lp1; ii < rows; ii++)
                            {
                                work[ii] += e[j] * matrix[j * rows + ii];
                            }
                        }

                        for (j = lp1; j < cols; j++)
                        {
                            double ww = -e[j] / e[lp1];
                            for (int ii = lp1; ii < rows; ii++)
                            {
                                matrix[j * rows + ii] += ww * work[ii];
                            }
                        }
                    }

                    if (computeVectors)
                    {
                        // place the transformation in v for subsequent back multiplication.
                        for (i = lp1; i < cols; i++)
                        {
                            v[l * cols + i] = e[i];
                        }
                    }
                }
            }

            //   set up the final bidiagonal matrix or order m.
            m = System.Math.Min(cols, rows + 1);
            nctp1 = nct + 1;
            nrtp1 = nrt + 1;
            if (nct < cols)
            {
                s[nctp1 - 1] = matrix[(nctp1 - 1) * rows + (nctp1 - 1)];
            }
            if (rows < m)
            {
                s[m - 1] = 0.0;
            }
            if (nrtp1 < m)
            {
                e[nrtp1 - 1] = matrix[(m - 1) * rows + (nrtp1 - 1)];
            }
            e[m - 1] = 0.0;

            //   if required, generate u.
            if (computeVectors)
            {
                for (j = nctp1 - 1; j < ncu; j++)
                {
                    for (i = 0; i < rows; i++)
                    {
                        u[j * rows + i] = 0.0;
                    }
                    u[j * rows + j] = 1.0;
                }

                for (l = nct - 1; l >= 0; l--)
                {
                    if (s[l] != 0.0)
                    {
                        for (j = l + 1; j < ncu; j++)
                        {
                            t = -Ddot(u, rows, l, j, l) / u[l * rows + l];
                            for (int ii = l; ii < rows; ii++)
                            {
                                u[j * rows + ii] += t * u[l * rows + ii];
                            }
                        }
                        DscalColumn(u, rows, l, l, -1.0);
                        u[l * rows + l] = 1.0 + u[l * rows + l];
                        for (i = 0; i < l; i++)
                        {
                            u[l * rows + i] = 0.0;
                        }
                    }
                    else
                    {
                        for (i = 0; i < rows; i++)
                        {
                            u[l * rows + i] = 0.0;
                        }
                        u[l * rows + l] = 1.0;
                    }
                }
            }

            //   if it is required, generate v.
            if (computeVectors)
            {
                for (l = cols - 1; l >= 0; l--)
                {
                    lp1 = l + 1;
                    if (l < nrt)
                    {
                        if (e[l] != 0.0)
                        {
                            for (j = lp1; j < cols; j++)
                            {
                                t = -Ddot(v, cols, l, j, lp1) / v[l * cols + lp1];
                                for (int ii = l; ii < cols; ii++)
                                {
                                    v[j * cols + ii] += t * v[l * cols + ii];
                                }
                            }
                        }
                    }
                    for (i = 0; i < cols; i++)
                    {
                        v[l * cols + i] = 0.0;
                    }
                    v[l * cols + l] = 1.0;
                }
            }

            //   transform s and e so that they are  double .
            for (i = 0; i < m; i++)
            {
                if (s[i] != 0.0)
                {
                    t = s[i];
                    r = s[i] / t;
                    s[i] = t;
                    if (i < m - 1)
                    {
                        e[i] = e[i] / r;
                    }
                    if (computeVectors)
                    {
                        DscalColumn(u, rows, i, 0, r);
                    }
                }
                //   ...exit
                if (i == m - 1)
                {
                    break;
                }
                if (e[i] != 0.0)
                {
                    t = e[i];
                    r = t / e[i];
                    e[i] = t;
                    s[i + 1] = s[i + 1] * r;
                    if (computeVectors)
                    {
                        DscalColumn(v, cols, i + 1, 0, r);
                    }
                }
            }

            //   main iteration loop for the singular values.
            mn = m;
            iter = 0;

            while (m > 0)
            { // quit if all the singular values have been found.
                // if too many iterations have been performed, set
                //      flag and return.
                if (iter >= MAXITER)
                {
                    info = m;
                    //   ......exit
                    break;
                }

                //      this section of the program inspects for
                //      negligible elements in the s and e arrays.  on
                //      completion the variables kase and l are set as follows.

                //         kase = 1     if s[m] and e[l-1] are negligible and l < m
                //         kase = 2     if s[l] is negligible and l < m
                //         kase = 3     if e[l-1] is negligible, l < m, and
                //                      s[l, ..., s[m] are not negligible (qr step).
                //         kase = 4     if e[m-1] is negligible (convergence).

                for (l = m - 2; l >= 0; l--)
                {
                    test = System.Math.Abs(s[l]) + System.Math.Abs(s[l + 1]);
                    ztest = test + System.Math.Abs(e[l]);
                    if (ztest == test)
                    {
                        e[l] = 0.0;
                        break;
                    }
                }

                if (l == m - 2)
                {
                    kase = 4;
                }
                else
                {
                    for (ls = m - 1; ls > l; ls--)
                    {
                        test = 0.0;
                        if (ls != m - 1)
                        {
                            test = test + System.Math.Abs(e[ls]);
                        }
                        if (ls != l + 1)
                        {
                            test = test + System.Math.Abs(e[ls - 1]);
                        }
                        ztest = test + System.Math.Abs(s[ls]);
                        if (ztest == test)
                        {
                            s[ls] = 0.0;
                            break;
                        }
                    }

                    if (ls == l)
                    {
                        kase = 3;
                    }
                    else if (ls == m - 1)
                    {
                        kase = 1;
                    }
                    else
                    {
                        kase = 2;
                        l = ls;
                    }
                }

                l = l + 1;

                //      perform the task indicated by kase.
                switch (kase)
                {
                    // deflate negligible s[m].
                    case 1:
                        f = e[m - 2];
                        e[m - 2] = 0.0;
                        for (k = m - 2; k >= 0; k--)
                        {
                            t1 = s[k];
                            Drotg(ref t1, ref f, ref cs, ref sn);
                            s[k] = t1;
                            if (k != l)
                            {
                                f = -sn * e[k - 1];
                                e[k - 1] = cs * e[k - 1];
                            }
                            if (computeVectors)
                            {
                                Drot(v, cols, k, m - 1, cs, sn);
                            }
                        }
                        break;

                    // split at negligible s[l].
                    case 2:
                        f = e[l - 1];
                        e[l - 1] = 0.0;
                        for (k = l; k < m; k++)
                        {
                            t1 = s[k];
                            Drotg(ref t1, ref f, ref cs, ref sn);
                            s[k] = t1;
                            f = -sn * e[k];
                            e[k] = cs * e[k];
                            if (computeVectors)
                            {
                                Drot(u, rows, k, l - 1, cs, sn);
                            }
                        }
                        break;

                    // perform one qr step.
                    case 3:
                        // calculate the shift.
                        scale = 0.0;
                        scale = System.Math.Max(scale, System.Math.Abs(s[m - 1]));
                        scale = System.Math.Max(scale, System.Math.Abs(s[m - 2]));
                        scale = System.Math.Max(scale, System.Math.Abs(e[m - 2]));
                        scale = System.Math.Max(scale, System.Math.Abs(s[l]));
                        scale = System.Math.Max(scale, System.Math.Abs(e[l]));
                        sm = s[m - 1] / scale;
                        smm1 = s[m - 2] / scale;
                        emm1 = e[m - 2] / scale;
                        sl = s[l] / scale;
                        el = e[l] / scale;
                        b = ((smm1 + sm) * (smm1 - sm) + emm1 * emm1) / 2.0;
                        c = (sm * emm1) * (sm * emm1);
                        shift = 0.0;
                        if (b != 0.0 || c != 0.0)
                        {
                            shift = System.Math.Sqrt(b * b + c);
                            if (b < 0.0)
                            {
                                shift = -shift;
                            }
                            shift = c / (b + shift);
                        }
                        f = (sl + sm) * (sl - sm) + shift;
                        g = sl * el;

                        // chase zeros.
                        for (k = l; k < m - 1; k++)
                        {
                            Drotg(ref f, ref g, ref cs, ref sn);
                            if (k != l) e[k - 1] = f;
                            f = cs * s[k] + sn * e[k];
                            e[k] = cs * e[k] - sn * s[k];
                            g = sn * s[k + 1];
                            s[k + 1] = cs * s[k + 1];
                            if (computeVectors)
                            {
                                Drot(v, cols, k, k + 1, cs, sn);
                            }
                            Drotg(ref f, ref g, ref cs, ref sn);
                            s[k] = f;
                            f = cs * e[k] + sn * s[k + 1];
                            s[k + 1] = -sn * e[k] + cs * s[k + 1];
                            g = sn * e[k + 1];
                            e[k + 1] = cs * e[k + 1];
                            if (computeVectors && k < rows)
                            {
                                Drot(u, rows, k, k + 1, cs, sn);
                            }
                        }
                        e[m - 2] = f;
                        iter = iter + 1;
                        break;

                    // convergence.
                    case 4:
                        // make the singular value  positive
                        if (s[l] < 0.0)
                        {
                            s[l] = -s[l];
                            if (computeVectors)
                            {
                                DscalColumn(v, cols, l, 0, -1.0);
                            }
                        }

                        // order the singular value.
                        while (l != mn - 1)
                        {
                            if (s[l] >= s[l + 1])
                            {
                                break;
                            }
                            t = s[l];
                            s[l] = s[l + 1];
                            s[l + 1] = t;
                            if (computeVectors && l < cols)
                            {
                                Dswap(v, cols, l, l + 1);
                            }
                            if (computeVectors && l < rows)
                            {
                                Dswap(u, rows, l, l + 1);
                            }
                            l = l + 1;
                        }
                        iter = 0;
                        m = m - 1;
                        break;
                }
            }
            if (computeVectors)
            {
                double[] temp = new double[cols * cols];
                System.Buffer.BlockCopy(v, 0, temp, 0, v.Length * Constants.SizeOfDouble);
                for (i = 0; i < cols; i++)
                {
                    for (j = 0; j < cols; j++)
                    {
                        v[i * cols + j] = temp[j * cols + i];
                    }
                }
            }
            return info;
        }

        #endregion Public Methods
        #region Private Methods

        private static double Dsign(double z1, double z2)
        {
            return System.Math.Abs(z1) * (z2 / System.Math.Abs(z2));
        }

        private static void Dswap(double[] A, int rows, int col1, int col2)
        {
            // swap columns Col1,Col2
            double z;
            for (int i = 0; i < rows; i++)
            {
                z = A[col1 * rows + i];
                A[col1 * rows + i] = A[col2 * rows + i];
                A[col2 * rows + i] = z;
            }
        }

        private static void DscalColumn(double[] A, int rows, int col, int start, double z)
        {
            // A part of column Col of matrix A from row Start to end multiply by z
            for (int i = start; i < rows; i++)
            {
                A[col * rows + i] = A[col * rows + i] * z;
            }
        }

        private static void DscalVector(double[] A, int start, double z)
        {
            // A part of vector  A from Start to end multiply by z
            for (int i = start; i < A.Length; i++)
            {
                A[i] = A[i] * z;
            }
        }

        private static void Drotg(ref double da, ref double db, ref double c, ref double s)
        {
            //     construct givens plane rotation.
            //     jack dongarra, linpack, 3/11/78.
            double roe, scale, r, z, Absda, Absdb, sda, sdb;

            roe = db;
            Absda = System.Math.Abs(da);
            Absdb = System.Math.Abs(db);
            if (Absda > Absdb) roe = da;
            scale = Absda + Absdb;
            if (scale == 0.0)
            {
                c = 1.0;
                s = 0.0;
                r = 0.0;
                z = 0.0;
            }
            else
            {
                sda = da / scale;
                sdb = db / scale;
                r = scale * System.Math.Sqrt(sda * sda + sdb * sdb);
                if (roe < 0.0) r = -r;
                c = da / r;
                s = db / r;
                z = 1.0;
                if (Absda > Absdb) z = s;
                if (Absdb >= Absda && c != 0.0) z = 1.0 / c;
            }
            da = r;
            db = z;
            return;
        }

        private static double Dnrm2Column(double[] A, int rows, int col, int start)
        {
            double s = 0;
            for (int i = start; i < rows; i++)
            {
                s += (A[col * rows + i] * A[col * rows + i]);
            }
            return System.Math.Sqrt(s);
        }

        private static double Dnrm2Vector(double[] A, int start)
        {
            double s = 0;
            for (int i = start; i < A.Length; i++)
            {
                s += (A[i] * A[i]);
            }
            return System.Math.Sqrt(s);
        }

        private static double Ddot(double[] A, int rows, int col1, int col2, int start)
        {
            double z = 0.0;
            for (int i = start; i < rows; i++)
            {
                z += A[col2 * rows + i] * A[col1 * rows + i];
            }
            return z;
        }

        private static void Drot(double[] A, int rows, int col1, int col2, double c, double s)
        {
            double z;
            for (int i = 0; i < rows; i++)
            {
                z = c * A[col1 * rows + i] + s * A[col2 * rows + i];
                A[col2 * rows + i] = c * A[col2 * rows + i] - s * A[col1 * rows + i];
                A[col1 * rows + i] = z;
            }
        }
        #endregion Private Methods
    }
}